COMMENT ⊗   VALID 00007 PAGES 
RECORD PAGE   DESCRIPTION
 00001 00001
 00003 00002	ENTRY INVERSION,BOUNDED,DECOMPOSE,SOLVE,IMPROVE
 00004 00003	INTERNAL SIMPLE PROCEDURE DECOMPOSE(INTEGER NSAFE REAL ARRAY A,LU)
 00008 00004	INTERNAL SIMPLE PROCEDURE SOLVE(INTEGER NSAFE REAL ARRAY LU,B,X)
 00010 00005	INTERNAL SIMPLE PROCEDURE IMPROVE(INTEGER NSAFE REAL ARRAY A,LU,B,XREFERENCE REAL DIGITS)
 00012 00006	INTERNAL SIMPLE PROCEDURE INVERSION(REAL ARRAY RES,MAT)
 00013 00007	INTERNAL BOOLEAN SIMPLE PROCEDURE BOUNDED(REAL X,YSAFE REAL ARRAY PINTEGER N)
 00016 ENDMK
⊗;
ENTRY INVERSION,BOUNDED,DECOMPOSE,SOLVE,IMPROVE;

BEGIN "SUPPORT"
SAFE INTEGER ARRAY PS[1:50];
SAFE REAL ARRAY R[1:50],DX[1:50];


SIMPLE PROCEDURE SINGULAR(INTEGER WHY);
	COMMENT PRINTS ERROR MESSAGES FOR DECOMPOSE AND IMPROVE;
BEGIN
	IF (WHY=0) THEN
	  USERERR(0,1, "MATRIX WITH ZERO ROW IN DECOMPOSE." );
	IF (WHY=1) THEN
	  USERERR(0,1, "SINGULAR MATRIX IN DECOMPOSE. SOLVE WILL DIVIDE BY ZERO." );
	IF (WHY=2) THEN
	  USERERR(0,1, "NO CONVERGENCE IN IMPROVE. MATRIX IS NEARLY SINGULAR." );
END ;
INTERNAL SIMPLE PROCEDURE DECOMPOSE(INTEGER N;SAFE REAL ARRAY A,LU);
		COMMENT A,LU[1:N,1:N];
		COMMENT USES GLOBAL INTEGER ARRAY PS;
		COMMENT COMPUTES TRIANGULAR MATRICES L AND U AND PERMUTATION
			MATRIX P SO THAT LU=PA. STORES L-I AND U IN LU.
			ARRAY PS CONTAINS PERMUTED ROW INDICES;
		COMMENT DECOMPOSE(N,A,A) OVERWRITES A WITH LU;
	BEGIN
		LABEL ENDKLOOP;
		INTEGER I,J,K,PIVOTINDEX;
		REAL NORMROW,PIVOT,SIZE,BIGGEST,MULT;
	SIMPLE PROCEDURE ILOOP(INTEGER UL;REFERENCE REAL R1,R2);
	START_CODE	LABEL LP,EU;
		MOVE 1,-1('17);
		MOVE 2,-2('17);
		MOVE 3,-3('17);
		SUB 3,K;
		JUMPLE 3,EU;
	LP:	AOJ 1,;
		AOJ 2,;
		MOVN 4,MULT;
		FMPR 4,(1);
		FADRM 4,(2);
		SOJG 3,LP;
	EU:	END;

		COMMENT INITIALIZE PS,LU AND R;
		FOR I←1 STEP 1 UNTIL N DO
			BEGIN
			PS[I]←I;
			NORMROW←0;
			FOR J←1 STEP 1 UNTIL N DO
			BEGIN
				LU[I,J]←A[I,J];
				IF (NORMROW<ABS(LU[I,J])) THEN NORMROW←ABS(LU[I,J]);
			END;
			IF (NORMROW≠0) THEN R[I]←1/NORMROW
			   ELSE BEGIN R[I]←0; SINGULAR(0); END;
			END;
			COMMENT GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING;
			FOR K←1 STEP 1 UNTIL N-1 DO
			BEGIN
				BIGGEST←0;
				FOR I←K STEP 1 UNTIL N DO
				BEGIN
				  SIZE←ABS(LU[PS[I],K])*R[PS[I]];
				  IF (BIGGEST<SIZE) THEN
					BEGIN BIGGEST←SIZE; PIVOTINDEX←I; END;
				END;
				IF BIGGEST=0 THEN
				  BEGIN SINGULAR(1); GO TO ENDKLOOP; END;
				IF PIVOTINDEX≠K THEN
				  BEGIN
				    J←PS[K];PS[K]←PS[PIVOTINDEX];PS[PIVOTINDEX]←J;
				  END;
				PIVOT←LU[PS[K],K];
				FOR I←K+1 STEP 1 UNTIL N DO
				BEGIN
				  LU[PS[I],K]←MULT←(LU[PS[I],K]/PIVOT);
				  IF MULT ≠0 THEN
COMMENT			  FOR J←K+1 STEP 1 UNTIL N DO
					LU[PS[I],J]←LU[PS[I],J]-MULT*LU[PS[K],J];
			ILOOP(N,LU[PS[I],K],LU[PS[K],K]);
					COMMENT INNER LOOP. ONLY COLUMN SUBSCRIPT
					  VARIES. USE MACHINE CODE IF NECESSARY
					  FOR EFFICIENCY;
				END;
ENDKLOOP:
END;
IF (LU[PS[N],N]=0) THEN SINGULAR(1);
END ;


INTERNAL SIMPLE PROCEDURE SOLVE(INTEGER N;SAFE REAL ARRAY LU,B,X);
	COMMENT LU[1:N,1:N],B,X[1:N];
	COMMENT USES GLOBAL SAFE INTEGER ARRAY PS;
	COMMENT SOLVES AX=B USING LU FROM DECOMPOSE;
BEGIN
	INTEGER I,J;
	REAL DOT;
	SIMPLE PROCEDURE ILOOP(INTEGER LL,UL;REFERENCE REAL R1,R2);
	START_CODE	LABEL LP,EU;
		MOVE 1,-1('17);
		MOVE 2,-2('17);
		MOVE 3,-3('17);
		SUB 3,-4('17);
		SETZ 4,;
		JUMPL 3,EU;
	LP:	MOVE 5,(1);
		FMPR 5,(2);
		FADR 4,5;
		AOJ 1,;
		AOJ 2,;
		SOJGE 3,LP;
	EU:	MOVEM 4,DOT;
	END;

	FOR I←1 STEP 1 UNTIL N DO
	BEGIN
COMMENT		DOT←0;
COMMENT		FOR J←1 STEP 1 UNTIL I-1 DO
		  DOT←DOT+LU[PS[I],J]*X[J];
		ILOOP(1,I-1,LU[PS[I],1],X[1]);
		X[I]←B[PS[I]]-DOT;
	END;
	FOR I←N STEP -1 UNTIL 1 DO
	BEGIN
COMMENT		DOT←0;
COMMENT		FOR J←I+1 STEP 1 UNTIL N DO
		  DOT←DOT+LU[PS[I],J]*X[J];
		ILOOP(I+1,N,LU[PS[I],I+1],X[I+1]);
		X[I]←(X[I]-DOT)/LU[PS[I],I];
	END;
	COMMENT AS IN DECOMPOSE, THE INNER LOOPS INVOLVE ONLY THE COLUMN
		SUBSCRIPT OF LU AND MAY BE MACHINE CODED FOR EFFICIENCY;
END;



INTERNAL SIMPLE PROCEDURE IMPROVE(INTEGER N;SAFE REAL ARRAY A,LU,B,X;REFERENCE REAL DIGITS);
	COMMENT A,LU[1:N,1:N],B,X[1:N];
	COMMENT A IS THE ORIGINAL MATRIX, LU IS FROM DECOMPOSE,B IS THE
		RIGHT HAND SIDE,X IS THE SOLUTION FROM SOLVE. IMPROVES
		X TO MACHINE ACCURACY AND SETS DIGITS TO THE NUMBER
		OF DIGITS OF X WHICH DO NOT CHANGE;
COMMENT MACHINE DEPENDENT QUANTITIES INDICATED BY 0-0;
BEGIN
		LABEL CONVERGED;
	INTEGER ITER, ITMAX,I;
	REAL RT,T,NORMX,NORMDX,EPS;
	EXTERNAL REAL SIMPLE PROCEDURE ADPFOR(INTEGER N;REAL ARRAY A;INTEGER I;REAL ARRAY X;REAL EXTRATERM);
	EPS←1.0@-8;
	ITMAX←16;
	NORMX←0;
	FOR I←1 STEP 1 UNTIL N DO
		IF (NORMX<ABS(X[I])) THEN NORMX←ABS(X[I]);
	IF NORMX=0 THEN GO TO CONVERGED;
	FOR ITER ←1 STEP 1 UNTIL ITMAX DO
	BEGIN
		FOR I←1 STEP 1 UNTIL N DO
			R[I]←ADPFOR(N,A,I,X,B[I]);
		SOLVE(N,LU,R,DX);
		NORMDX←0;
		FOR I←1 STEP 1 UNTIL N DO
		BEGIN
			T←X[I];
			X[I]←X[I]+DX[I];
			IF (NORMDX<ABS(X[I]-T)) THEN NORMDX←ABS(X[I]-T);
		END;
		IF (NORMDX≤EPS*NORMX) THEN GO TO CONVERGED;
	END ;
	COMMENT ITERATION DID NOT CONVERGE;
	SINGULAR(2);
CONVERGED:
END;
INTERNAL SIMPLE PROCEDURE INVERSION(REAL ARRAY RES,MAT);
BEGIN "INVERT"
SAFE OWN REAL ARRAY LU[1:4,1:4],B,X[1:4];
INTEGER I,J;
REAL DIGITS;
DECOMPOSE(4,MAT,LU);
FOR J←1 STEP 1 UNTIL 4 DO BEGIN
	FOR I← 1 STEP 1 UNTIL 4 DO B[I]←IF I=J THEN 1.0 ELSE 0.0;
	SOLVE(4,LU,B,X);
	IMPROVE(4,MAT,LU,B,X,DIGITS);
	FOR I←1 STEP 1 UNTIL 4 DO RES[I,J]←X[I] END;
END "INVERT";
INTERNAL BOOLEAN SIMPLE PROCEDURE BOUNDED(REAL X,Y;SAFE REAL ARRAY P;INTEGER N);
BEGIN	INTEGER I,C,N2,II;
	REAL M1,M2,RV,RH;
	LABEL NEXTL,ADD4;
	C←0;
	FOR I←N STEP 2 UNTIL 2*(N-1) DO
	BEGIN	RV←(P[(I+1) MOD N]-Y)*(Y-P[(I+3) MOD N]);
		IF RV<0.0 THEN GO TO NEXTL;
		IF RV> 0.0 THEN 
		BEGIN	RH←(P[I MOD N]-X)*(X-P[(I+2) MOD N]);
			IF RH < 0.0 THEN
			BEGIN	C←IF X<P[I MOD N] THEN C+1 ELSE C;
				GO TO NEXTL
			END;
			IF RH>0.0 THEN
			BEGIN	M1←(P[(I+1) MOD N]-P[(I+3) MOD N])/
				(P[I MOD N]-P[(I+2) MOD N]);
				M2←(IF P[(I+1) MOD N]>P[(I+3) MOD N]  THEN
				(P[(I+1) MOD N]-Y)/(P[I MOD N]-X) ELSE
				(P[(I+3) MOD N]-Y)/(P[(I+2) MOD N]-X)) ;
				IF M1=M2 THEN RETURN (TRUE);
				C←IF M1>M2 THEN C+1 ELSE C;
				GO TO NEXTL
			END;
			IF P[I MOD N]=P[(I+2) MOD N] THEN RETURN (TRUE);
			IF P[I MOD N]=X THEN
			BEGIN	C←IF X<P[(I+2) MOD N] THEN C+1 ELSE C
			END ELSE C←IF X<P[I MOD N] THEN C+1 ELSE C;
			GO TO NEXTL
		END;
		IF P[(I+1) MOD N]=P[(I+3) MOD N] THEN BEGIN IF (P[I MOD N]-X)*(X-P[(I+2) MOD N])≥0.0 
		THEN RETURN (TRUE) ELSE GO TO NEXTL END;
		II← IF Y=P[(I+1) MOD N] THEN I ELSE I+2;
		IF X>P[II MOD N] THEN GO TO ADD4;
		IF X<P[II MOD N] THEN
		BEGIN	C←IF (P[(II+3) MOD N]-P[(II+1) MOD N])*
			(P[(II+1) MOD N]-P[(II-1) MOD N])>0.0 THEN C+1 ELSE C;
			GO TO ADD4
		END;
		RETURN (TRUE);
	ADD4:	I←I+2;
	NEXTL:	END;
	RETURN(IF C  MOD 2 =1 THEN TRUE ELSE FALSE)
END;
END "SUPPORT"